      subroutine wclbes(zz,eta1,zlmin,nl,fc,gc,fcp,gcp,sig,kfn,mode1,
     1                  ifail,ipr)
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c                                                                      c
c  complex coulomb wavefunction program using steed's method           c
c  original title : coulcc                                             c
c                                                                      c
c  a. r. barnett           manchester  march   1981                    c
c  modified i.j. thompson  daresbury, sept. 1983 for complex functions c
c                                                                      c
c  the form (not the substance) of this program has been modified      c
c   by k.s. kolbig (cern)    december 1987                             c
c                                                                      c
c  original program  rcwfn       in    cpc  8 (1974) 377-395           c
c                 +  rcwff       in    cpc 11 (1976) 141-142           c
c                 +  coulfg      in    cpc 27 (1982) 147-166           c
c  description of real algorithm in    cpc 21 (1981) 297-314           c
c  description of complex algorithm    jcp 64 (1986) 490-509           c
c  this version written up       in    cpc 36 (1985) 363-372           c
c                                                                      c
c  wclbes returns f,g,g',g',sig for complex eta, zz, and zlmin,        c
c   for nl integer-spaced lambda values zlmin to zlmin+nl inclusive,   c
c   thus giving  complex-energy solutions to the coulomb schrodinger   c
c   equation,to the klein-gordon equation and to suitable forms of     c
c   the dirac equation ,also spherical & cylindrical bessel equations  c
c                                                                      c
c  if abs(mode1)                                                       c
c            = 1  get f,g,f',g'   for integer-spaced lambda values     c
c            = 2      f,g      unused arrays must be dimensioned in    c
c            = 3      f,  f'          call to at least length (0:1)    c
c            = 4      f                                                c
c            = 11 get f,h+,f',h+' ) if kfn=0, h+ = g + i.f        )    c
c            = 12     f,h+        )       >0, h+ = j + i.y = h(1) ) in c
c            = 21 get f,h-,f',h-' ) if kfn=0, h- = g - i.f        ) gc c
c            = 22     f,h-        )       >0, h- = j - i.y = h(2) )    c
c                                                                      c
c     if mode1 < 0 then the values returned are scaled by an exponen-  c
c                  tial factor (dependent only on zz) to bring nearer  c
c                  unity the functions for large abs(zz), small eta ,  c
c                  and abs(zl) < abs(zz).                              c
c        define scale = (  0        if mode1 > 0                       c
c                       (  imag(zz) if mode1 < 0  &  kfn < 3           c
c                       (  real(zz) if mode1 < 0  &  kfn = 3           c
c        then fc = exp(-abs(scale)) * ( f, j, j, or i)                 c
c         and gc = exp(-abs(scale)) * ( g, y, or y )                   c
c               or exp(scale)       * ( h+, h(1), or k)                c
c               or exp(-scale)      * ( h- or h(2) )                   c
c                                                                      c
c  if  kfn  =  0,-1  complex coulomb functions are returned   f & g    c
c           =  1   spherical bessel      "      "     "       j & y    c
c           =  2 cylindrical bessel      "      "     "       j & y    c
c           =  3 modified cyl. bessel    "      "     "       i & k    c
c                                                                      c
c          and where coulomb phase shifts put in sig if kfn=0 (not -1) c
c                                                                      c
c  the use of mode and kfn is independent                              c
c    (except that for kfn=3,  h(1) & h(2) are not given)               c
c                                                                      c
c  with negative orders lambda, wclbes can still be used but with      c
c    reduced accuracy as cf1 becomes unstable. the user is thus        c
c    strongly advised to use reflection formulae based on              c
c    h+-(zl,,) = h+-(-zl-1,,) * exp +-i(sig(zl)-sig(-zl-1)-(zl+1/2)pi) c
c                                                                      c
c  precision:  results to within 2-3 decimals of 'machine accuracy',   c
c              except in the following cases:                          c
c              (1) if cf1a fails because x too small or eta too large  c
c               the f solution  is less accurate if it decreases with  c
c               decreasing lambda (e.g. for lambda.le.-1 & eta.ne.0)   c
c              (2) if eta is large (e.g. >> 50) and x inside the       c
c                turning point, then progressively less accuracy       c
c              (3) if zlmin is around sqrt(accur) distance from an     c
c               integral order and abs(x) << 1, then errors present.   c
c               rerr traces the main roundoff errors.                  c
c                                                                      c
c   wclbes is coded for real*8 on ibm or equivalent  accur >= 10**-14  c
c          with a section of doubled real*16 for less roundoff errors. c
c          (if no doubled precision available, increase jmax to eg 100)c
c   use implicit complex*32 & real*16 on vs compiler accur >= 10**-32  c
c   for single precision cdc (48 bits) reassign                        c
c        double precision=real etc.                                    c
c                                                                      c
c   ipr    on input   = 0 : no printing of error messages              c
c                    ne 0 : print error messages on file 6             c
c   ifail  in output = -2 : argument out of range                      c
c                    = -1 : one of the continued fractions failed,     c
c                           or arithmetic check before final recursion c
c                    =  0 : all calculations satisfactory              c
c                    ge 0 : results available for orders up to & at    c
c                             position nl-ifail in the output arrays.  c
c                    = -3 : values at zlmin not found as over/underflowc
c                    = -4 : roundoff errors make results meaningless   c
c                                                                      c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c                                                                      c
c     machine dependent constants :                                    c
c                                                                      c
c     accu     target bound on relative error (except near 0 crossings)c
c               (accur should be at least 100 * acc8)                  c
c     acc8    smallest number with 1+acc8 .ne.1 in real*8  arithmetic  c
c     acc16    smallest number with 1+acc16.ne.1 in real*16 arithmetic c
c     fpmax    magnitude of largest floating point number * acc8       c
c     fpmin    magnitude of smallest floating point number / acc8      c
c                                                                      c
c     parameters determining region of calculations :                  c
c                                                                      c
c        r20      estimate of (2f0 iterations)/(cf2 iterations)        c
c        asym     minimum x/(eta**2+l) for cf1a to converge easily     c
c        xnear    minimum abs(x) for cf2 to converge accurately        c
c        limit    maximum no. iterations for cf1, cf2, and 1f1 series  c
c        jmax     size of work arrays for pade accelerations           c
c        ndrop    number of successive decrements to define instabilityc
c                                                                      c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c     c309r1 = cf2,   c309r2 = f11,    c309r3 = f20,
c     c309r4 = cf1r,  c309r5 = cf1c,   c309r6 = cf1a,
c     c309r7 = rcf,   c309r8 = ctidy.
c
      implicit complex*16 (a-h,o-z)
c
      common /coulc2/ nfp,n11,npq1,npq2,n20,kas1,kas2
      double precision zero,one,two,half
      double precision r20,asym,xnear
      double precision fpmax,fpmin,fplmax,fplmin
      double precision accu,acc8,acc16
      double precision hpi,tlog
      double precision err,rerr,absc,accur,acct,acch,accb,c309r4
      double precision paccq,eps,off,scale,sf,sfsh,ta,rk,omega,absx
 
      logical lpr,etane0,ifcp,rlel,donem,unstab,zlneg,axial,nocf2,npint
 
      parameter(zero = 0, one = 1, two = 2, half = one/two)
      parameter(ci = (0,1), cih = half*ci)
      parameter(r20 = 3, asym = 3, xnear = half)
      parameter(limit = 20000, ndrop = 5, jmax = 50)
 
      parameter(fpmax = 3d290, fpmin = 2d-300)
      parameter(accu = 1d-14, acc8 = 2d-16, acc16 = 2d-16)
 
      parameter(hpi  = 1.57079 63267 94896 619d0)
      parameter(tlog = 0.69314 71805 59945 309d0)
 
      dimension fc(0:*),gc(0:*),fcp(0:*),gcp(0:*),sig(0:*)
      dimension xrcf(jmax,4)
      integer npq(2),kas(2)
      equivalence (npq(1),npq1),(npq(2),npq2)
      equivalence (kas(1),kas1),(kas(2),kas2)
 
      nintc(w)=nint(dreal(w))
      absc(w)=abs(dreal(w))+abs(dimag(w))
      npint(w,accb)=absc(nintc(w)-w) .lt. accb .and. dreal(w) .lt. half
c
      mode=mod(abs(mode1),10)
      ifcp=mod(mode,2) .eq. 1
      lpr=ipr .ne. 0
      ifail=-2
      n11=0
      nfp=0
      kas(1)=0
      kas(2)=0
      npq(1)=0
      npq(2)=0
      n20=0
 
      accur=max(accu,50*acc8)
      acct=half*accur
      acch=sqrt(accur)
      accb=sqrt(acch)
      rerr=acct
      fplmax=log(fpmax)
      fplmin=log(fpmin)
c
      cik=one
      if(kfn .ge. 3) cik=ci*sign(one,fpmin-dimag(zz))
      x=zz*cik
      eta=eta1
      if(kfn .gt. 0) eta=zero
      etane0=absc(eta) .gt. acc8
      etai=eta*ci
      dell=zero
      if(kfn .ge. 2) dell=half
      zm1=zlmin-dell
      scale=zero
      if(mode1 .lt. 0) scale=dimag(x)
c
      m1=1
      l1=m1+nl
      rlel=abs(dimag(eta))+abs(dimag(zm1)) .lt. acc8
      absx=abs(x)
      axial=rlel .and. abs(dimag(x)) .lt. acc8*absx
      if(mode .le. 2 .and. absx .lt. fpmin) go to 310
      xi=one/x
      xlog=log(x)
 
c       log with cut along the negative real axis, see also omega
 
      id=1
      donem=.false.
      unstab=.false.
      lf=m1
      ifail=-1
   10 zlm=zm1+(lf-m1)
      zll=zm1+(l1-m1)
c
c ***       zll  is final lambda value, or 0.5 smaller for j,y bessels
c
      z11=zll
      if(id .lt. 0) z11=zlm
      p11=ci*sign(one,acc8-dimag(eta))
      last=l1
c
c ***       find phase shifts and gamow factor at lambda = zll
c
      pk=zll+one
      aa=pk-etai
      ab=pk+etai
      bb=two*pk
      zlneg=npint(bb,accb)
      clgaa=wlogam(aa)
      clgab=clgaa
      if(etane0 .and. .not.rlel) clgab=wlogam(ab)
      if(etane0 .and. rlel) clgab=dconjg(clgaa)
      sigma=(clgaa-clgab)*cih
      if(kfn .eq. 0) sig(l1)=sigma
      if(.not.zlneg) cll=zll*tlog-hpi*eta-wlogam(bb)+(clgaa+clgab)*half
      theta=x-eta*(xlog+tlog)-zll*hpi+sigma
c
      ta=(dimag(aa)**2+dimag(ab)**2+abs(dreal(aa))+abs(dreal(ab)))*half
      if(id .gt. 0 .and. absx .lt. ta*asym .and. .not.zlneg) go to 20
c
c ***         use cf1 instead of cf1a, if predicted to converge faster,
c                 (otherwise using cf1a as it treats negative lambda &
c                  recurrence-unstable cases properly)
c
      rk=sign(one,dreal(x)+acc8)
      p=theta
      if(rk .lt. 0) p=-x+eta*(log(-x)+tlog)-zll*hpi-sigma
      xrcf(1,1)=pk
      f=rk*c309r6(x*rk,eta*rk,zll,p,acct,jmax,nfp,fest,err,fpmax,xrcf,
     1                                      xrcf(1,3),xrcf(1,4))
      fesl=log(fest)+abs(dimag(x))
      nfp=-nfp
      if(nfp .lt. 0 .or. unstab .and. err .lt. accb) go to 40
      if(.not.zlneg .or. unstab .and. err .gt. accb) go to 20
      if(lpr) write(6,1060) '-l',err
      if(err.gt.accb) go to 280
      go to 40
c
c ***    evaluate cf1  =  f   =  f'(zll,eta,x)/f(zll,eta,x)
c
   20 if(axial) then
       f=c309r4(x,eta,zll,acc8,sf,rk,etane0,limit,err,nfp,
     1          fpmin,fpmax,lpr)
       fcl=sf
       tpk1=rk
      else
       f=c309r5(x,eta,zll,acc8,fcl,tpk1,etane0,limit,err,nfp,
     1          fpmin,fpmax,lpr)
      end if
      if(err .gt. one) then
       ifail=-1
       go to 290
      end if
c
c ***  make a simple check for cf1 being badly unstable:
c
      if(id .ge. 0) then
       unstab=dreal((one-eta*xi)*ci*dimag(theta)/f) .gt. zero
     1  .and. .not.axial .and. abs(dimag(theta)) .gt. -log(acc8)*half
     2  .and. absc(eta)+absc(zll) .lt. absc(x)
       if(unstab) then
        id=-1
        lf=l1
        l1=m1
        rerr=acct
        go to 10
       end if
      end if
c
c *** compare accumulated phase fcl with asymptotic phase for g(k+1) :
c     to determine estimate of f(zll) (with correct sign) to start recur
c
      w=x*x*(half/tpk1+one/tpk1**2)+eta*(eta-two*x)/tpk1
      fesl=(zll+one)*xlog+cll-w-log(fcl)
   40 fesl=fesl-abs(scale)
      rk=max(dreal(fesl),fplmin*half)
      fesl=dcmplx(min(rk,fplmax*half),dimag(fesl))
      fest=exp(fesl)
c
      rerr=max(rerr,err,acc8*abs(dreal(theta)))
c
      fcl=fest
      fpl=fcl*f
      if(ifcp) fcp(l1)=fpl
      fc(l1)=fcl
c
c *** downward recurrence to lambda = zlm. array gc,if present,stores rl
c
      i=max(-id,0)
      zl=zll+i
      mono=0
      off=abs(fcl)
      ta=absc(sigma)
 
c
c     corresponds to   do 70 l = l1-id,lf,-id
c
      l=l1-id
      lc70=(l1-lf)/id
   70 if(lc70 .le. 0) go to 80
      if(etane0) then
       if(rlel) then
        dsig=atan2(dreal(eta),dreal(zl))
        rl=sqrt(dreal(zl)**2+dreal(eta)**2)
       else
        aa=zl-etai
        bb=zl+etai
        if(absc(aa) .lt. acch .or. absc(bb) .lt. acch) goto 50
        dsig=(log(aa)-log(bb))*cih
        rl=aa*exp(ci*dsig)
       end if
       if(absc(sigma) .lt. ta*half) then
 
c               re-calculate sigma because of accumulating roundoffs:
 
        sl=(wlogam(zl+i-etai)-wlogam(zl+i+etai))*cih
        rl=(zl-etai)*exp(ci*id*(sigma-sl))
        sigma=sl
        ta=zero
       else
        sigma=sigma-dsig*id
       end if
       ta=max(ta,absc(sigma))
       sl=eta+zl*zl*xi
       pl=zero
       if(absc(zl) .gt. acch) pl=(sl*sl-rl*rl)/zl
       fcl1=(fcl*sl+id*zl*fpl)/rl
       sf=abs(fcl1)
       if(sf .gt. fpmax) go to 350
       fpl=(fpl*sl+id*pl*fcl)/rl
       if(mode .le. 1) gcp(l+id)=pl*id
      else
 
c                      eta = 0, including bessels.  nb rl==sl
 
       rl=zl*xi
       fcl1=fcl*rl+fpl*id
       sf=abs(fcl1)
       if(sf .gt. fpmax) go to 350
       fpl=(fcl1*rl-fcl)*id
      end if
      mono=mono+1
      if(sf .ge. off) mono=0
      fcl=fcl1
      off=sf
      fc(l)=fcl
      if(ifcp) fcp(l)=fpl
      if(kfn .eq. 0) sig(l)=sigma
      if(mode .le. 2) gc(l+id)=rl
      zl=zl-id
      if(mono .lt. ndrop .or. axial .or.
     1             dreal(zlm)*id .gt. -ndrop .and. .not.etane0) then
       l=l-id
       lc70=lc70-1
       go to 70
      end if
      unstab=.true.
c
c ***    take action if cannot or should not recur below this zl:
 
   50 zlm=zl
      lf=l
      if(id .lt. 0) go to 380
      if(.not.unstab) lf=l+1
      if(l+mono .lt. l1-2 .or. id .lt. 0 .or. .not.unstab) go to 80
 
c             otherwise, all l values (for stability) should be done
c                        in the reverse direction:
 
      id=-1
      lf=l1
      l1=m1
      rerr=acct
      go to 10
 
   80 if(fcl .eq. zero) fcl=acc8
      f=fpl/fcl
c
c *** check, if second time around, that the 'f' values agree!
c
      if(id .gt. 0) first=f
      if(donem) rerr=max(rerr,absc(f-first)/absc(f))
      if(donem) go to 90
c
      nocf2=.false.
      thetam=x-eta*(xlog+tlog)-zlm*hpi+sigma
c
c *** on left x-plane, determine omega by requiring cut on -x axis
c     on right x-plane, choose omega (using estimate based on thetam)
c       so h(omega) is smaller and recurs upwards accurately.
c     (x-plane boundary is shifted to give cf2(lh) a chance to converge)
c
      omega=sign(one,dimag(x)+acc8)
      if(dreal(x) .ge. xnear) omega=sign(one,dimag(thetam)+acc8)
      sfsh=exp(omega*scale-abs(scale))
      off=exp(min(two*max(abs(dimag(x)),abs(dimag(thetam)),
     1                         abs(dimag(zlm))*3),fplmax))
      eps=max(acc8,acct*half/off)
c
c ***    try first estimated omega, then its opposite,
c        to find the h(omega) linearly independent of f
c        i.e. maximise  cf1-cf2 = 1/(f h(omega)) , to minimise h(omega)
c
   90 do 100 l = 1,2
      lh=1
      if(omega .lt. zero) lh=2
      pm=ci*omega
      etap=eta*pm
      if(donem) go to 130
      pq1=zero
      paccq=one
      kase=0
c
c ***            check for small x, i.e. whether to avoid cf2 :
c
      if(mode .ge. 3 .and. absx .lt. one ) go to 190
      if(mode .lt. 3 .and. (nocf2 .or. absx .lt. xnear .and.
     1   absc(eta)*absx .lt. 5 .and. absc(zlm) .lt. 4)) then
       kase=5
       go to 120
      end if
c
c ***  evaluate   cf2 : pq1 = p + i.omega.q  at lambda = zlm
c
      pq1=c309r1(x,eta,zlm,pm,eps,limit,err,npq(lh),acc8,acch,
     1             lpr,accur,dell)
      err=err*max(one,absc(pq1)/max(absc(f-pq1),acc8))
c
c *** check if impossible to get f-pq accurately because of cancellation
c                original guess for omega (based on thetam) was wrong
c                use kase 5 or 6 if necessary if re(x) < xnear
      if(err .lt. acch) go to 110
      nocf2=dreal(x) .lt. xnear .and. abs(dimag(x)) .lt. -log(acc8)
  100 omega=-omega
      if(unstab) go to 360
      if(lpr .and. dreal(x) .lt. -xnear) write(6,1060) '-x',err
  110 rerr=max(rerr,err)
c
c ***  establish case of calculation required for irregular solution
c
  120 if(kase .ge. 5) go to 130
      if(dreal(x) .gt. xnear) then
 
c          estimate errors if kase 2 or 3 were to be used:
 
       paccq=eps*off*absc(pq1)/max(abs(dimag(pq1)),acc8)
      end if
      if(paccq .lt. accur) then
       kase=2
       if(axial) kase=3
      else
       kase=1
       if(npq(1)*r20 .lt. jmax) kase=4
c             i.e. change to kase=4 if the 2f0 predicted to converge
      end if
  130 go to (190,140,150,170,190,190), abs(kase)
  140 if(.not.donem)
c
c ***  evaluate   cf2 : pq2 = p - i.omega.q  at lambda = zlm   (kase 2)
c
     1  pq2=c309r1(x,eta,zlm,-pm,eps,limit,err,npq(3-lh),acc8,acch,
     2             lpr,accur,dell)
      p=(pq2+pq1)*half
      q=(pq2-pq1)*half*pm
      go to 160
  150 p=dreal(pq1)
      q=dimag(pq1)
c
c ***   with kase = 3 on the real axes, p and q are real & pq2 = pq1*
c
      pq2=dconjg(pq1)
c
c *** solve for fcm = f at lambda = zlm,then find norm factor w=fcm/fcl
c
  160 w=(pq1-f)*(pq2-f)
      sf=exp(-abs(scale))
      fcm=sqrt(q/w)*sf
 
c                  any sqrt given here is corrected by
c                  using sign for fcm nearest to phase of fcl
 
      if(dreal(fcm/fcl) .lt. zero) fcm=-fcm
      gam=(f-p)/q
      ta=absc(gam+pm)
      paccq=eps*max(ta,one/ta)
      hcl=fcm*(gam+pm)*(sfsh/(sf*sf))
 
c                            consider a kase = 1 calculation
 
      if(paccq .gt. accur .and. kase .gt. 0) then
       f11v=c309r2(x,eta,z11,p11,acct,limit,0,err,n11,fpmax,acc8,acc16)
       if(err .lt. paccq) go to 200
      end if
      rerr=max(rerr,paccq)
      go to 230
c
c *** arrive here if kase = 4
c     to evaluate the exponentially decreasing h(lh) directly.
c
  170 if(donem) go to 180
      aa=etap-zlm
      bb=etap+zlm+one
      f20v=c309r3(aa,bb,-half*pm*xi,acct,jmax,err,fpmax,n20,xrcf)
      if(n20 .le. 0) go to 190
      rerr=max(rerr,err)
      hcl=fpmin
      if(abs(dreal(pm*thetam)+omega*scale) .gt. fplmax) go to 330
  180 hcl=f20v*exp(pm*thetam+omega*scale)
      fcm=sfsh/((f-pq1)*hcl)
      go to 230
c
c *** arrive here if kase=1   (or if 2f0 tried mistakenly & failed)
c
c           for small values of x, calculate f(x,sl) directly from 1f1
c               using dreal*16 arithmetic if possible.
c           where z11 = zll if id>0, or = zlm if id<0
c
  190 f11v=c309r2(x,eta,z11,p11,acct,limit,0,err,n11,fpmax,acc8,acc16)
  200 if(n11 .lt. 0) then
 
c                               f11 failed from bb = negative integer
 
       if(lpr) write(6,1060) '-l',one
       ifail=-1
       go to 290
      end if
 
c                      consider a kase 2 or 3 calculation :
 
      if(err .gt. paccq .and. paccq .lt. accb) then
       kase=-2
       if(axial) kase=-3
       go to 130
      end if
      rerr=max(rerr,err)
      if(err .gt. fpmax) go to 370
      if(id .lt. 0) cll=z11*tlog-hpi*eta-wlogam(bb)
     1                        +wlogam(z11+one+p11*eta)-p11*sigma
      ek=(z11+one)*xlog-p11*x+cll-abs(scale)
      if(id .gt. 0) ek=ek-fesl+log(fcl)
      if(dreal(ek) .gt. fplmax) go to 350
      if(dreal(ek) .lt. fplmin) go to 340
      fcm=f11v*exp(ek)
      if(kase .ge. 5) then
       if(absc(zlm+zlm-nintc(zlm+zlm)) .lt. acch) kase=6
c
c ***  for abs(x) < xnear, then cf2 may not converge accurately, so
c ***      use an expansion for irregular soln from origin :
c
       sl=zlm
       zlneg=dreal(zlm) .lt. -one+accb
       if(kase .eq. 5 .or. zlneg) sl=-zlm-one
       pk=sl+one
       aa=pk-etap
       ab=pk+etap
       bb=two*pk
       clgaa=wlogam(aa)
       clgab=clgaa
       if(etane0) clgab=wlogam(ab)
       clgbb=wlogam(bb)
       if(kase .eq. 6 .and. .not.zlneg) then
        if(npint(aa,accur)) clgaa=clgab-two*pm*sigma
        if(npint(ab,accur)) clgab=clgaa+two*pm*sigma
       end if
       cll=sl*tlog-hpi*eta-clgbb+(clgaa+clgab)*half
       dsig=(clgaa-clgab)*pm*half
       if(kase .eq. 6) p11=-pm
       ek=pk*xlog-p11*x+cll-abs(scale)
       sf=exp(-abs(scale))
       chi=zero
       if(.not.(kase .eq. 5 .or. zlneg)) go to 210
c
c ***  use  g(l)  =  (cos(chi) * f(l) - f(-l-1)) /  sin(chi)
c
c      where chi = sig(l) - sig(-l-1) - (2l+1)*pi/2
c
       chi=sigma-dsig-(zlm-sl)*hpi
       f11v=c309r2(x,eta,sl,p11,acct,limit,0,err,npq(1),
     1             fpmax,acc8,acc16)
       rerr=max(rerr,err)
       if(kase .eq. 6) go to 210
       fesl=f11v*exp(ek)
       fcl1=exp(pm*chi)*fcm
       hcl=fcl1-fesl
       rerr=max(rerr,acct*max(absc(fcl1),absc(fesl))/absc(hcl))
       hcl=hcl/sin(chi)*(sfsh/(sf*sf))
       go to 220
c
c *** use the logarithmic expansion for the irregular solution (kase 6)
c        for the case that bb is integral so sin(chi) would be zero.
c
  210  rl=bb-one
       n=nintc(rl)
       zlog=xlog+tlog-pm*hpi
       chi=chi+pm*thetam+omega*scale+ab*zlog
       aa=one-aa
       if(npint(aa,accur)) then
        hcl=zero
       else
        if(id .gt. 0 .and. .not.zlneg) f11v=fcm*exp(-ek)
        hcl=exp(chi-clgbb-wlogam(aa))*(-1)**(n+1)*(f11v*zlog+
     1   c309r2(x,eta,sl,-pm,acct,limit,2,err,npq(2),fpmax,acc8,acc16))
        rerr=max(rerr,err)
       end if
       if(n .gt. 0) then
        ek=chi+wlogam(rl)-clgab-rl*zlog
        df=c309r2(x,eta,-sl-one,-pm,zero,n,0,err,l,fpmax,acc8,acc16)
        hcl=hcl+exp(ek)*df
       end if
       rerr=max(rerr,two*abs(bb-nintc(bb)))
  220  pq1=f-sfsh/(fcm*hcl)
      else
       if(mode .le. 2) hcl=sfsh/((f-pq1)*fcm)
       kase=1
      end if
c
c ***  now have absolute normalisations for coulomb functions
c          fcm & hcl  at lambda = zlm
c      so determine linear transformations for functions required :
c
  230 ih=abs(mode1)/10
      if(kfn .eq. 3) ih=(3-dimag(cik))/2+half
      p11=one
      if(ih .eq. 1) p11=ci
      if(ih .eq. 2) p11=-ci
      df=-pm
      if(ih .ge. 1) df=-pm+p11
      if(absc(df) .lt. acch) df=zero
c
c *** normalisations for spherical or cylindrical bessel functions
c
      if(kfn .le. 0) then
       alfa=zero
       beta=one
      else if(kfn .eq. 1) then
       alfa=xi
       beta=xi
      else
       alfa=xi*half
       beta=sqrt(xi/hpi)
       if(dreal(beta) .lt. zero) beta=-beta
      end if
      aa=one
      if(kfn .gt. 0) aa=-p11*beta
 
c                calculate rescaling factors for i & k output
 
      if(kfn .ge. 3) then
       p=exp((zlm+dell)*hpi*cik)
       aa=beta*hpi*p
       beta=beta/p
       q=cik*id
      end if
 
c                  calculate rescaling factors for gc output
 
      if(ih .eq. 0) then
       ta=abs(scale)+dimag(pm)*scale
       rk=zero
       if(ta .lt. fplmax) rk=exp(-ta)
      else
       ta=abs(scale)+dimag(p11)*scale
       if(absc(df) .gt. acch .and. ta .gt. fplmax) go to 320
       if(absc(df) .gt. acch) df=df*exp(ta)
       sf=two*(lh-ih)*scale
       rk=zero
       if(sf .gt. fplmax) go to 320
       if(sf .gt. fplmin) rk=exp(sf)
      end if
      kas((3-id)/2)=kase
      w=fcm/fcl
      if(log(absc(w))+log(absc(fc(lf))) .lt. fplmin) go to 340
      if(mode .ge. 3) go to 240
      if(lpr .and. absc(f-pq1) .lt. acch*absc(f))
     1                             write(6,1020) lh,zlm+dell
      hpl=hcl*pq1
      if(absc(hpl) .lt. fpmin .or. absc(hcl) .lt. fpmin) go to 330
c
c *** idward recurrence from hcl,hpl(lf) (stored gc(l) is rl if reqd)
c *** renormalise fc,fcp at each lambda
c ***    zl   = zlm - min(id,0) here
c
  240 do 270 l = lf,l1,id
      fcl=w*fc(l)
      if(absc(fcl) .lt. fpmin) go to 340
      if(ifcp) fpl=w*fcp(l)
      fc(l)=beta*fcl
      if(ifcp) fcp(l)=beta*(fpl-alfa*fcl)*cik
      fc(l)=c309r8(fc(l),accur)
      if(ifcp) fcp(l)=c309r8(fcp(l),accur)
      if(mode .ge. 3) go to 260
      if(l .eq. lf) go to 250
      zl=zl+id
      zid=zl*id
      rl=gc(l)
      if(etane0) then
       sl=eta+zl*zl*xi
       if(mode .eq. 1) then
        pl=gcp(l)
       else
        pl=zero
        if(absc(zl) .gt. acch) pl=(sl*sl-rl*rl)/zid
       end if
       hcl1=(sl*hcl-zid*hpl)/rl
       hpl=(sl*hpl-pl*hcl)/rl
      else
       hcl1=rl*hcl-hpl*id
       hpl=(hcl-rl*hcl1)*id
      end if
      hcl=hcl1
      if(absc(hcl) .gt. fpmax) go to 320
  250 gc(l)=aa*(rk*hcl+df*fcl)
      if(mode .eq. 1) gcp(l)=(aa*(rk*hpl+df*fpl)-alfa*gc(l))*cik
      gc(l)=c309r8(gc(l),accur)
      if(mode .eq. 1) gcp(l)=c309r8(gcp(l),accur)
      if(kfn .ge. 3) aa=aa*q
  260 if(kfn .ge. 3) beta=-beta*q
  270 last=min(last,(l1-l)*id)
      go to 280
c
c *** come here after all soft errors to determine how many l values ok
c
  310 if(lpr) write(6,1000) zz
      go to 999
  320 if(lpr) write(6,1010) zl+dell,'ir',hcl,'>',fpmax
      go to 280
  330 if(lpr) write(6,1010) zl+dell,'ir',hcl,'<',fpmin
      go to 280
  340 if(lpr) write(6,1010) zl+dell,'  ',fcl,'<',fpmin
      go to 280
  350 if(lpr) write(6,1010) zl+dell,'  ',fcl,'>',fpmax
      go to 280
  360 if(lpr) write(6,1030) zl+dell
      go to 280
  370 if(lpr) write(6,1040) z11,i
      ifail=-1
      go to 290
  380 if(lpr) write(6,1050) zlmin,zlm,zlm+one,zlmin+nl
      ifail=-1
      go to 290
  280 if(id .gt. 0 .or. last .eq. 0) ifail=last
      if(id .lt. 0 .and. last .ne. 0) ifail=-3
c
c *** come here after all errors for this l range (zlm,zll)
c
c *** so on first block, 'f' started decreasing monotonically,
c                        or hit bound states for low zl.
c     thus redo m1 to lf-1 in reverse direction, i.e. do
c      cf1a at zlmin & cf2 at zlm (midway between zlmin & zlmax)
c
  290 if(id .gt. 0 .and. lf .ne. m1) then
       id=-1
       if(.not.unstab) lf=lf-1
       donem=unstab
       lf=min(lf,l1)
       l1=m1
       go to 10
      end if
      if(ifail .lt. 0) go to 999
      if(lpr .and. rerr .gt. accb) write(6,1070) rerr
      if(rerr .gt. 0.1d0) ifail=-4
  999 do 998 l = 1,nl+1
      fc(l-1)=fc(l)
      gc(l-1)=gc(l)
      fcp(l-1)=fcp(l)
      gcp(l-1)=gcp(l)
  998 sig(l-1)=sig(l)
      return
c
 1000 format(1x,'***** cern c309 wclbes ... ',
     1 'cannot calculate irregular solutions for x =',
     2 1p,2d10.2,' abs(x) too small')
 1010 format(1x,'***** cern c309 wclbes ... ',
     1 'at zl =',2f8.3,' ',a2,'regular solution (',1p,2e10.1,') ',
     2 a1,e10.1)
 1020 format(1x,'***** cern c309 wclbes ... ',
     1 'warning: linear independence between ''f'' and ''h(',i1,
     2 ')'' is lost at zl = ',2f7.2/1x,'*****',22x,'(eg. coulomb ',
     3 'eigenstate or cf1 unstable)')
 1030 format(1x,'***** cern c309 wclbes ... ',
     2 '(eta & l)/x too large for cf1a, and cf1 unstable at l = ',2f8.2)
 1040 format(1x,'***** cern c309 wclbes ... ',
     1 'overflow in 1f1 series at zl = ',2f8.3,' at term',i5)
 1050 format(1x,'***** cern c309 wclbes ... ',
     1 'both bound-state poles and f-instabilities occur or multiple',
     2 ' instabilities present'/
     3   1x,'*****',22x,'try calling twice, 1st for zl from',2f8.3,
     4 ' to',2f8.3,' (incl)'/1x,'*****',41x,'2nd for zl from',2f8.3,
     5 ' to',2f8.3)
 1060 format(1x,'***** cern c309 wclbes ... ',
     1 'warning: as ''',a2,''' reflection rules not used ',
     2 'errors can be up to',1pd12.2)
 1070 format(1x,'***** cern c309 wclbes ... ',
     1 'warning: overall roundoff error approximately',1pe11.1)
      end
cdeck  id>, wdigam.
      function wdigam(u)
c*ns  complex*16 wdigam,w,u,v,h,r,p
c*ns  double precision c(6),x,y,a,pi,zero,one,half,c1,d
      complex*16 wdigam,  u,v,h,r,p
      double precision c(6),x,y,a,pi,zero,one,half,c1
 
      data pi /3.14159 26535 89793 24d0/
      data zero /0d0/, one /1d0/, half /0.5d0/, c1 /15d0/
 
      data c(1) / 8.33333 33333 33333 33d-2/
      data c(2) /-8.33333 33333 33333 33d-3/
      data c(3) / 3.96825 39682 53968 25d-3/
      data c(4) /-4.16666 66666 66666 67d-3/
      data c(5) / 7.57575 75757 57575 76d-3/
      data c(6) /-2.10927 96092 79609 28d-2/
 
c*ul    9 x=dreal(u)
      x=dreal(u)
      a=abs(x)
      if(dimag(u) .eq. zero .and. -a .eq. int(x)) then
       wdigam=zero
       write(6,100) x
       return
      end if
      if(x .lt. zero) u=-u
      v=u
      h=zero
      if(a .lt. c1) then
       h=one/v
       do 1 i = 1,14-int(a)
       v=v+one
    1  h=h+one/v
       v=v+one
      end if
      r=one/v**2
      p=r*c(1)
      do 2 i = 6,1,-1
    2 p=r*(c(i)+p)
      h=log(v)-half/v-p-h
      if(x .lt. zero) then
       v=pi*u
       x=dreal(v)
       y=tanh(dimag(v))
       a=sin(x)
       x=cos(x)
       h=h+one/u+pi*dcmplx(x,-a*y)/dcmplx(a,x*y)
      end if
      wdigam=h
 
  100 format(7x,'***** cern c307 wdigam ... argument equals ',
     1          'non-positive integer = ',e10.1)
      end
      complex function cdigam(z)
      complex z
      complex*16 wdigam,w,ww
      double precision d
 
      sround(d) = sngl(d+(d-dble(sngl(d))))
c*ia  d
      data d/0/
 
      w  = z
      ww = wdigam(w)
      cdigam = cmplx(sround(dreal(ww)),sround(dimag(ww)))
 
      end
cdeck  id>, wlogam.
      function wlogam(v)
c*ns  complex*16 wlogam,w,u,v,h,p,r
c*ns  double precision c(10),pi,t,x,y,ya,a,ur,ui,c1,c2,d
      complex*16 wlogam,  u,v,h,p,r
      double precision c(10),pi,t,x,y,ya,a,ur,ui,c1,c2
      double precision zero,one,seven,half
 
      data zero /0d0/, one /1d0/, seven /7d0/, half /0.5d0/
 
      data pi /3.14159 26535 89793 24d+0/
      data c1 /9.18938 53320 46727 42d-1/
      data c2 /1.14472 98858 49400 17d+0/
 
      data c( 1) / 8.33333 33333 33333 33d-2/
      data c( 2) /-2.77777 77777 77777 78d-3/
      data c( 3) / 7.93650 79365 07936 51d-4/
      data c( 4) /-5.95238 09523 80952 38d-4/
      data c( 5) / 8.41750 84175 08417 51d-4/
      data c( 6) /-1.91752 69175 26917 53d-3/
      data c( 7) / 6.41025 64102 56410 26d-3/
      data c( 8) /-2.95506 53594 77124 18d-2/
      data c( 9) / 1.79644 37236 88305 73d-1/
      data c(10) /-1.39243 22169 05901 12d+0/
 
c*ul9 x=dreal(v)
      x=dreal(v)
      y=dimag(v)
      if(y .eq. zero .and. -abs(x) .eq. int(x)) then
       wlogam=zero
       write(6,100) x
       return
      end if
      ya=abs(y)
      u=dcmplx(x,ya)
      if(x .lt. zero) u=one-u
      h=zero
      ur=dreal(u)
      if(ur .lt. seven) then
       ui=dimag(u)
       a=atan2(ui,ur)
       h=u
       do 1 i = 1,6-int(ur)
       ur=ur+one
       u=dcmplx(ur,ui)
       h=h*u
    1  a=a+atan2(ui,ur)
       h=dcmplx(half*log(dreal(h)**2+dimag(h)**2),a)
       u=u+one
      end if
      r=one/u**2
      p=r*c(10)
      do 2 i = 9,2,-1
    2 p=r*(c(i)+p)
      h=c1+(u-half)*log(u)-u+(c(1)+p)/u-h
      if(x .lt. zero) then
       ur=int(x)-one
       ui=pi*(x-ur)
       x=pi*ya
       t=exp(-x-x)
       a=sin(ui)
       t=x+half*log(t*a**2+(half*(one-t))**2)
       a=atan2(cos(ui)*tanh(x),a)-ur*pi
       h=c2-dcmplx(t,a)-h
      end if
      if(y .lt. zero) h=dconjg(h)
      wlogam=h
 
  100 format(7x,'***** cern c333 wlogam ... argument equals ',
     1          'non-positive integer = ',e10.1)
      end
      complex function clogam(z)
      complex z
      complex*16 wlogam,w,ww
      double precision d
 
      sround(d) = sngl(d+(d-dble(sngl(d))))
c*ia  d
      data d/0/
 
      w  = z
      ww = wlogam(w)
      clogam = cmplx(sround(dreal(ww)),sround(dimag(ww)))
 
      end
cdeck  id>, r1dp.
      function c309r1(x,eta,zl,pm,eps,limit,err,npq,acc8,acch,
     1             lpr,accur,dell)
c
c                                    (omega)        (omega)
c *** evaluate  cf2  = p + pm.q  =  h   (eta,x)' / h   (eta,x)
c                                    zl             zl
c     where pm = omega.i
c
      implicit complex*16(a-h,o-z)
      logical lpr
      double precision eps,err,acc8,acch,accur,ta,rk
      double precision absc,half
 
      parameter(half = 1d0/2d0)
 
      absc(w)=abs(dreal(w))+abs(dimag(w))
 
      ta=limit+limit
      etap=eta*pm
      xi=1/x
      wi=etap+etap
      rk=0
      pq=(1-eta*xi)*pm
      aa=-(eta*eta+zl*zl+zl)+etap
      bb=2*(x-eta+pm)
      rl=xi*pm
      if(absc(bb) .lt. acch) then
       rl=rl*aa/(aa+rk+wi)
       pq=pq+rl*(bb+pm+pm)
       aa=aa+2*(rk+1+wi)
       bb=bb+4*pm
       rk=rk+4
      end if
      dd=1/bb
      dl=aa*dd*rl
   10 pq=pq+dl
      rk=rk+2
      aa=aa+rk+wi
      bb=bb+pm+pm
      dd=1/(aa*dd+bb)
      dl=dl*(bb*dd-1)
      err=absc(dl)/absc(pq)
      if(err .ge. max(eps,acc8*rk*half) .and. rk .le. ta) go to 10
c
      npq=half*rk
      c309r1=pq+dl
      if(lpr .and. npq .ge. limit-1 .and. err .gt. accur)
     1 write(6,1000) int(dimag(pm)),npq,err,zl+dell
      return
 1000 format(1x,'***** cern c309 wclbes ... ',
     2 'cf2(',i2,') not converged fully in ',i7,' iterations'/1x,27x,
     3 'error in irregular solution =',1p,d11.2,' at zl = ',2f8.3)
      end
cdeck  id>, r2dp.
      function c309r2(x,eta,zl,p,eps,limit,kind,err,nits,
     1                fpmax,acc8,acc16)
c
c *** evaluate the hypergeometric function
c                                        i
c            f (aa;bb; z) = sum  (aa)   z / ( (bb)  i! )
c           1 1              i       i            i
c
c     to accuracy eps with at most limit terms.
c  if kind = 0 : using extended precision but real arithmetic only,
c            1 : using normal precision in complex arithmetic,
c   or       2 : using normal complex arithmetic, but with wdigam factor
c
c  where aa, bb, and z are defined below
c
      implicit double precision(a-h,o-z)
      complex*16 x,eta,zl,p,aa,bb,z,c309r2,wdigam
      complex*16 dd,g,f,ai,bi,t
      double precision ar,br,gr,gi,dr,di,tr,ti,ur,ui,fi,fi1,den
 
      parameter(tbbb = 3d0/2d0)
 
      absc(aa)=abs(dreal(aa))+abs(dimag(aa))
      nintc(aa)=nint(dreal(aa))
       c309r2=0
      aa=zl+1-eta*p
      bb=2*zl+2
      z=2*p*x
      if(dreal(bb) .le. 0 .and. abs(bb-nintc(bb)) .lt.
     1 sqrt(sqrt(acc8))**3 .and. dreal(bb)+limit .ge. tbbb) then
       nits=-1
       return
      end if
      if(limit .le. 0) then
       c309r2=0
       err=0
       nits=1
       return
      end if
      ta=1
      rk=1
      if(kind .le. 0 .and. absc(z)*absc(aa) .gt. absc(bb)*1.0) then
       dr=1
       di=0
       gr=1
       gi=0
       ar=dreal(aa)
       br=dreal(bb)
       fi=0
       do 20 i = 2,limit
       fi1=fi+1
       tr=br*fi1
       ti=dimag(bb)*fi1
       den=1/(tr*tr+ti*ti)
       ur=(ar*tr+dimag(aa)*ti)*den
       ui=(dimag(aa)*tr-ar*ti)*den
       tr=ur*gr-ui*gi
       ti=ur*gi+ui*gr
       gr=dreal(z)*tr-dimag(z)*ti
       gi=dreal(z)*ti+dimag(z)*tr
       dr=dr+gr
       di=di+gi
       err=abs(gr)+abs(gi)
       if(err .gt. fpmax) go to 60
       rk=abs(dr)+abs(di)
       ta=max(ta,rk)
       if(err .lt. rk*eps .or. i .ge. 4 .and. err .lt. acc16) go to 30
       fi=fi1
       ar=ar+1
   20  br=br+1
c
   30  c309r2=dr+(0,1)*di
       err=acc16*ta/rk
      else
c
c*    if real*16 arithmetic is not available, (or already using it!),
c*    then use kind > 0
c
       g=1
       f=1
       if(kind .ge. 2) f=wdigam(aa)-wdigam(bb)-wdigam(g)
       dd=f
       do 40 i = 2,limit
       ai=aa+(i-2)
       bi=bb+(i-2)
       r=i-1
       g=g*z*ai/(bi*r)
c
c                       multiply by (psi(a+r)-psi(b+r)-psi(1+r))
c
       if(kind .eq. 2) f=f+1/ai-1/bi-1/r
       t=g*f
       dd=dd+t
       err=absc(t)
       if(err .gt. fpmax) go to 60
       rk=absc(dd)
       ta=max(ta,rk)
       if(err .lt. rk*eps .or. err .lt. acc8 .and. i .ge. 4) go to 50
   40  continue
 
   50  err=acc8*ta/rk
       c309r2=dd
      end if
   60 nits=i
      return
      end
cdeck  id>, r3dp.
      function c309r3(aa,bb,z,eps,jmax,re,fpmax,n,x)
c
c     evaluate the hypergeometric function
c                                             i
c            f (aa,bb;;z) = sum  (aa)  (bb)  z / i!
c           2 0              i       i     i
c
c     to accuracy eps with at most jmax terms.
c
c     if the terms start diverging,
c     the corresponding continued fraction is found by rcf
c     & evaluated progressively by steed's method to obtain convergence.
c
      implicit complex*16(a-h,o-z)
      dimension x(jmax,4)
      logical finite
      double precision ep,eps,at,atl,absc,re,fpmax
 
      absc(w)=abs(dreal(w))+abs(dimag(w))
      nintc(w)=nint(dreal(w))
c
      re=0
      x(1,1)=1
      sum=x(1,1)
      atl=absc(x(1,1))
      f=sum
      d=1
      df=sum
      j=0
      ep=eps*(10*jmax)
      ma=-nintc(aa)
      mb=-nintc(bb)
      finite=abs(abs(dreal(aa))-ma) .lt. ep .and. abs(dimag(aa)) .lt. ep
     1  .or. abs(abs(dreal(bb))-mb) .lt. ep .and. abs(dimag(bb)) .lt. ep
      imax=jmax
      if(finite .and. ma .ge. 0) imax=min(ma+1,imax)
      if(finite .and. mb .ge. 0) imax=min(mb+1,imax)
      do 10 i = 2,imax
      x(i,1)=x(i-1,1)*z*(aa+i-2)*(bb+i-2)/(i-1)
      if(absc(x(i,1)) .gt. fpmax) then
       n=0
       c309r3=sum
       if(.not.finite) re=at/absc(sum)
       return
      end if
      at=absc(x(i,1))
      if(j .eq. 0) then
       sum=sum+x(i,1)
       if(at .lt. absc(sum)*eps) then
        n=i
        c309r3=sum
        if(.not.finite) re=at/absc(sum)
        return
       end if
      end if
      if(finite) go to 10
      if(j .gt. 0 .or. at .gt. atl .or. i .ge. jmax-2) j=j+1
      if(j .eq. 0) go to 10
      call c309r7(x(1,1),x(1,2),j,i,x(1,3),eps)
      if(i .lt. 0) then
       n=0
       c309r3=sum
       if(.not.finite) re=at/absc(sum)
       return
      end if
      do 50 k = max(j,2),i
      d=1/(d*x(k,2)+1)
      df=df*d-df
      f=f+df
      if(absc(df) .lt. absc(f)*eps .or.
     1   df .eq. 0 .and. f .eq. 0 .and. i .ge. 4) then
       n=k
       c309r3=f
       re=absc(df)/absc(f)
       return
      end if
   50 continue
      j=i
   10 atl=at
      if(.not.finite) i=-jmax
      n=i
      c309r3=sum
      if(.not.finite) re=at/absc(sum)
      return
      end
cdeck  id>, r4dp.
      function c309r4(x,eta,zl,eps,fcl,tpk1,etane0,limit,err,nfp,
     1              fpmin,fpmax,lpr)
c
c ***   evaluate cf1 = f'(zl,eta,x)/f(zl,eta,x)    (real)
c
      implicit double precision(a-h,o-z)
      logical lpr,etane0
 
      fcl=1
      xi=1/x
      pk=zl+1
      px=pk+limit
      ek=eta/pk
      f=ek+pk*xi
      if(abs(f) .lt. fpmin) f=fpmin
      d=0
      c=f
      small=sqrt(fpmin)
      rk2=1+ek*ek
c
c ***   begin cf1 loop on pk = k = lambda + 1
c
   10 pk1=pk+1
      tpk1=pk+pk1
      if(etane0) then
       ek=eta/pk
       rk2=1+ek*ek
       tk=tpk1*(xi+ek/pk1)
      else
       tk=tpk1*xi
      end if
      c=tk-rk2/c
      d=tk-rk2*d
      if(abs(c) .lt. fpmin) c=fpmin
      if(abs(d) .lt. fpmin) d=fpmin
      d=1/d
      df=d*c
      f=f*df
      fcl=fcl*d*tpk1*xi
      if(abs(fcl) .lt. small) fcl=fcl/small
      if(abs(fcl) .gt. fpmax) fcl=fcl*fpmin
      pk=pk1
      if(pk .le. px) then
       if(abs(df-1) .ge. eps) go to 10
       nfp=pk-zl-1
       err=eps*sqrt(real(nfp))
       c309r4=f
      else
       if(lpr) write (6,1000) limit,abs(x)
       err=2
      end if
      return
 1000 format(1x,'***** cern c309 wclbes ... cf1 (real) has failed ',
     1'to converge after',i10,' iterations as abs(x) =',f15.0)
      end
cdeck  id>, r5dp.
      function c309r5(x,eta,zl,eps,fcl,tpk1,etane0,limit,err,nfp,
     1              fpmin,fpmax,lpr)
c
c ***    evaluate cf1  =  f'(zl,eta,x)/f(zl,eta,x)  (complex)
c
      implicit complex*16(a-h,o-z)
      logical lpr,etane0
      double precision eps,err,fpmin,fpmax,absc,small,px
 
 
      absc(w)=abs(dreal(w))+abs(dimag(w))
 
      fcl=1
      xi=1/x
      pk=zl+1
      px=pk+limit
      ek=eta/pk
      f=ek+pk*xi
      if(absc(f) .lt. fpmin) f=fpmin
      d=0
      c=f
      small=sqrt(fpmin)
      rk2=1+ek*ek
c
c ***   begin cf1 loop on pk = k = lambda + 1
c
   10 pk1=pk+1
      tpk1=pk+pk1
      if(etane0) then
       ek=eta/pk
       rk2=1+ek*ek
       tk=tpk1*(xi+ek/pk1)
      else
       tk=tpk1*xi
      end if
      c=tk-rk2/c
      d=tk-rk2*d
      if(absc(c) .lt. fpmin) c=fpmin
      if(absc(d) .lt. fpmin) d=fpmin
      d=1/d
      df=d*c
      f=f*df
      fcl=fcl*d*tpk1*xi
      if(absc(fcl) .lt. small) fcl=fcl/small
      if(absc(fcl) .gt. fpmax) fcl=fcl*fpmin
      pk=pk1
      if(dreal(pk) .le. px) then
       if(absc(df-1) .ge. eps) go to 10
       nfp=pk-zl-1
       err=eps*sqrt(real(nfp))
       c309r5=f
      else
       if(lpr) write (6,1000) limit,abs(x)
       err=2
      end if
      return
 1000 format(1x,'***** cern c309 wclbes ... cf1 (complex) has failed ',
     1'to converge after',i10,' iterations as abs(x) =',f15.0)
      end
cdeck  id>, r6dp.
      function c309r6(rho,eta,xl,psi,eps,nmax,nused,fcl,re,fpmax,xx,g,c)
c
c     evaluate the asymptotic expansion for the
c            logarithmic derivative of the regular solution
c
c ***        cf1a  =  f   =  f'(xl,eta,rho)/f(xl,eta,rho)
c
c      that is valid for dreal(rho)>0, and best for rho >> eta**2, xl,
c      and is derived from the 2f0 expansions for h+ and h-
c      e.g. by froeberg (rev. mod. physics vol 27, p399 , 1955)
c      some lines of this subprogram are for convenience copied from
c           takemasa, tamura & wolter cpc 17 (1979) 351.
c
c     evaluate to accuracy eps with at most nmax terms.
c
c     if the terms start diverging,
c     the corresponding continued fraction is found by rcf
c     & evaluated progressively by steed's method to obtain convergence.
c
      implicit complex*16(a-h,o-z)
      dimension xx(2,nmax),g(nmax),c(nmax)
      double precision re,eps,t1,t2,t3,at,atl,absc,fpmax
      double precision hpi
 
      parameter(hpi  = 1.57079 63267 94896 619d0)
 
      absc(w)=abs(dreal(w))+abs(dimag(w))
c
      t1=sin(dreal(psi))
      t2=cos(dreal(psi))
      atl=tanh(dimag(psi))
 
c         give cos(psi)/cosh(im(psi)), which always has correct sign
 
      cosl=dcmplx(t2,-t1*atl)
      tanl=dcmplx(t1,t2*atl)/cosl
      re=0
      xll1=xl*xl+xl
      etasq=eta*eta
      sl1=1
      sl=sl1
      sc1=0
      sc=sc1
      tl1=sc
      tl=tl1
      tc1=1-eta/rho
      tc=tc1
      fcl=tl+sl*tanl
      g(1)=(tc+sc*tanl)/fcl
      glast=g(1)
      atl=absc(glast)
      f=glast
      d=1
      df=glast
      j=0
      do 10 n = 2,nmax
      t1=n-1
      t2=2*t1-1
      t3=t1*t1-t1
      denom=2*rho*t1
      c1=(eta*t2)/denom
      c2=(etasq+xll1-t3)/denom
      sl2=c1*sl1-c2*tl1
      tl2=c1*tl1+c2*sl1
      sc2=c1*sc1-c2*tc1-sl2/rho
      tc2=c1*tc1+c2*sc1-tl2/rho
      sl=sl+sl2
      tl=tl+tl2
      sc=sc+sc2
      tc=tc+tc2
      sl1=sl2
      tl1=tl2
      sc1=sc2
      tc1=tc2
      fcl=tl+sl*tanl
      if(absc(fcl) .gt. fpmax .or. absc(fcl) .lt. 1./fpmax) then
       c309r6=g(1)
       fcl=1
       re=1
       nused=0
       return
      end if
      gsum=(tc+sc*tanl)/fcl
      g(n)=gsum-glast
      glast=gsum
      at=absc(g(n))
      if(at .lt. absc(gsum)*eps) then
       fcl=fcl*cosl
       c309r6=gsum
       re=at/absc(gsum)
       nused=n
       return
      end if
      if(j .gt. 0 .or. at .gt. atl .or. n .ge. nmax-2) j=j+1
      if(j .eq. 0) go to 10
      call c309r7(g,c,j,n,xx,eps)
      if(n .lt. 0) then
       c309r6=g(1)
       fcl=1
       re=1
       nused=0
       return
      end if
      do 60 k = max(j,2),n
      d=1/(d*c(k)+1)
      df=df*d-df
      f=f+df
      if(absc(df) .lt. absc(f)*eps .or.
     1   df .eq. 0 .and. f .eq. 0 .and. n .ge. 4) then
       c309r6=f
       fcl=fcl*cosl
       re=absc(df)/absc(f)
       nused=k
       return
      end if
   60 continue
      j=n
   10 atl=at
      c309r6=f
      fcl=fcl*cosl
      re=absc(df)/absc(f)
      nused=-nmax
      return
      end
cdeck  id>, r7dp.
      subroutine c309r7(a,b,ibeg,inum,xx,eps)
c
c*******************************************************************
c
c  rcf converts polynomial a to the corresponding continued
c         fraction, in 'normal'  form with coefficients b
c         by the 'p algorithm' of patry & gupta
c
c   a(z) = a1/z + a2/z**3 + a3/z**5 + ... + an/z**(2n-1)
c
c   b(z) = b1/z+ b2/z+ b3/z+ .../(z+ bn/z)
c
c  data:
c   a     vector a(k), k=1,inum         input
c   b     vector b(k), k=ibeg,inum      output
c   ibeg  order of first coef. calc.    input
c   inum  order of a, even or odd       input
c   xx    auxiliary vector of length .ge. length of vector b
c         caller provides space for a,b,xx
c   note that neither of the first two terms a(1) a(2) should be zero
c          & the user can start the calculation with any value of
c          ibeg provided the c.f. coefs have been already
c          calculated up to inum = ibeg-1
c          & the method breaks down as soon as the absolute value
c          of a c.f. coef. is less than eps.    at the time of the
c          break up  inum has been replaced by minus times the number
c          of this coefficient.
c   algorithm: j. patry & s. gupta, eir-bericht 247, november 1973
c              eidg. institut fur reaktorforschung
c              wuerenlingen, switzerland
c   see also:  haenggi, roesel & trautmann,
c              j. comput. phys., v. 137, (1980) 242-258
c   note:      restart procedure modified by i.j.thompson
c
c*******************************************************************
c
      implicit complex*16(a-h,o-z)
      dimension a(100),b(100),xx(2,100)
      logical even
      double precision eps
 
      ibn=inum
      if(ibeg .gt. 4) go to 50
      if(ibeg .eq. 4) go to 20
      b(1)=a(1)
      if(ibn .ge. 2) b(2)=-a(2)/a(1)
      if(ibn .lt. 3) return
      x0=a(3)/a(2)
      xx(2,1)=b(2)
      xx(1,1)=-x0
      xx(1,2)=0
      b(3)=-x0-b(2)
      x0=-b(3)*a(2)
      m=3
      mp12=2
      even=.true.
      if(ibn .le. 3) return
   20 if(abs(b(3)) .lt. eps*abs(x0)) then
       inum=-m
       return
      end if
      m=4
   30 x1=a(m)
      m2m1=mp12
      mp12=m2m1+1
      if(even) mp12=m2m1
      do 40 k = 2,mp12
   40 x1=x1+a(m-k+1)*xx(1,k-1)
      b(m)=-x1/x0
      if(m .ge. ibn) return
   50 if(abs(b(m)) .lt. eps*abs(x0)) then
       inum=-m
       return
      end if
      do 60 k = m2m1,2,-1
   60 xx(2,k)=xx(1,k)+b(m)*xx(2,k-1)
      xx(2,1)=xx(1,1)+b(m)
      do 70 k = 1,m2m1
      x0=xx(2,k)
      xx(2,k)=xx(1,k)
   70 xx(1,k)=x0
      x0=x1
      xx(1,m2m1+1)=0
      m=m+1
      even=.not.even
      go to 30
      end
cdeck  id>, r8dp.
      function c309r8(z,acc)
      complex*16 c309r8,z
      double precision acc
      double precision x,y,ax,ay,a
 
      x=dreal(z)
      y=dimag(z)
      ax=abs(x)
      ay=abs(y)
      a=5*acc*(ax+ay)
      if(ax .lt. a) x=0
      if(ay .lt. a) y=0
      c309r8=dcmplx(x,y)
      return
      end
